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The basic models of helix-coil transitions in biomolecules are introduced. 
These include phenomenological, zipper (Bragg-Zimm) models of polypeptides, 
loop-entropy (Poland-Scheraga) and Hamiltonian (Peyrard-Bishop) models of ho- 
mogeneous DNA denaturation. The transfer integral approach to one-dimensional 
thermodynamics is presented in some detail, including the necessary extensions to 
deal with the singular integral equations arising in the case of on-site potentials 
with a flat top. The (non-) applicability of the theorems which prohibit phase 
transitions in one-dimensional systems is discussed. 

1 Introduction 

The purpose of these notes is to provide a brief introduction to concepts and 
methods employed in the description of thermodynamic phase transitions in 
model biomolecular systems. This is one of the areas where the LOCNET 
network has been active, with an emphasis on modelling the fundamental 
interactions which provide a basis for understanding both the cooperative 
behavior and the nonlinear dynamics (what I will call "Hamiltonian" mod- 
els). My primary aim is therefore to present a working introduction to this 
area, including some necessary details on methods and tools; I hope that 
young researchers who enter the field will find this material useful. 

On the other hand there is a long and distinguished tradition in the field, 
which has achieved remarkable progress, based on a phenomenological de- 
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scription of the statistical properties of helix formation and growth and the 
entropies associated with them. I have therefore chosen to include a section 
on fundamentals of these "zipper" and "loop" -based models Again, the 
hope is that researchers who are active in developing microscopic Hamilto- 
nian models will be guided by the succesful features of phenomenological, 
"zipper" and "loop" models. In addition, there is a simple, utilitarian rea- 
son for doing this: most of the experimental data presently available has 
been analyzed, directly or indirectly, in terms of such models. Therefore, 
anyone seriously interested in comparing theory and experiment should be 
familiar with them. 

The plan of this review is as follows: Section 2 deals with zipper and 
loop models of polypeptides and DNA, respectively. Section 3 describes 
the Hamiltonian approach to DNA denaturation, including a somewhat de- 
tailed introduction to the transfer integral method. An appendix discusses 
the (non-) applicability of theorems which prohibit phase transitions in one 
dimension to the models presented. 



2 Zippers and Loops 

2.1 Helix-Coil transitions in polypeptides 
2.1.1 Background 

Synthetic polypeptides, i.e. macromolecules consisting of identical amino- 
acid residues, are ideal for studying the transition from the alpha-helical 
to coil-like structure. Understanding of this transition is central to con- 
trolling the stability of secondary protein structure |Q. Residues in helical 
regions give rise to distinct experimental signatures (e.g. viscocity, optical 
rotation). At a given macromolecular size N - which can be controlled in 
synthetic polypeptides - one can measure the helix fraction as a function 
of temperature. Typically^, that fraction completes the transition from 1 
to over a fairly narrow temperature range - a few degrees K in the case 
of long chains. Chemists describe the process A{helix) < — > B{coil) as an 
equilibrium between the two species, 

K^^^ (1) 

CA 

where the helix fraction is given by 

e^^^ = ^. (2) 

CA + CB l + K 

The sign convention is as follows: I am looking at the conversion of helix 
(A) to coil (B); therefore AG^Gb-Ga^AH- TAS, and AH > 0, i.e. 
the helix is energetically favored. 



2 




270 280 290 300 
Temperature (K) 



Figure 1: Coil fraction vs. temperature for a polypeptide (poly-7-benzyl- 
L-glutamate) of controlled {N = 1500) length; The curve has been calcu- 
lated in the framework of the generalized zipper model - cf. Eq. |l^ below 
- (redrawn after Ref. Q ). 

The value O = 0.5 defines the midpoint of the transition, Tm- Assuming 
(although this is not exact, and sometimes not even a good approximation) 
that the enthalpy and entropy differences do not depend very much on 
temperature, leads to 



dQ 
dT 



e=o.5 



1 AH 



(3) 



The inverse of Eq. (g) measures the width of the transition (in degrees K). 
A sharp transition (of a few degrees K) has a high [van't Hoff] AH (of the 
order of 100 Kcal/mol), indicating that perhaps as many as 100 hydrogen 
bonds are cooperatively broken during the transition. 



2.1.2 "Zipper" model family: underlying concepts 

Helix initiation and helix growth are viewedH] as distinct processes: 
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• Growth: an existing helix may grow further at the nth site, or shrink. 
This is viewed as a forward and reverse reaction, with a rate ratio 
s = exp(— AG*/T), which reflects the difference in local free energies 
between the helix and coil states. If the ratio is greater than unity, 
the helix has a tendency to grow ("zip"). If it is less than unity, 
the helix will shrink ("unzip"). At temperatures near the transition, 
s « 1. The enthalpy difference AH* < corresponds to the energy 
of a single hydrogen bond formed in the process of helix growth. 

• Nucleation: in order to initiate a helix, 3 residues have to organize 
themselves. Again, viewing nucleation as a forward / reverse reaction, 
introduces a dimensionless a = exp{—AGinit/T). The large difference 
in the free energy comes mostly from the entropy loss associated with 
the organization of the 3-4 residues involved in the first turn of the 
helix. 

I now present an outline of theoretical models [Q, in order of increasing 
complexity: 



2.1.3 "0-th order" - The "all or nothing" (AON) model: 

Only two states are significant within this model. The pure coil, with 
relative statistical weight equal to unity; and the helix with N residues, 
with a relative weight as^ . Intermediate states are suppressed, presumably 
due to high rate barriers. This gives a helix fraction 

and a slope at midpoint 

'de\ N AH* 



"-^ I e=o.5 ^ ^ni 



(5) 



There is strong cooperativity. 



2.1.4 Further considerations: the zipper model 

The model allows a single connected helical region of any length n < N . 
The statistical weight (Boltzmann factor) is -according to the general con- 
siderations, cf. above-, cts", and the helix can commence at any of the first 
An = N — n + 1 positions. This gives a partition function 

N 

Z = l + J2 ^"'^■s" (6) 

n=l 
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and a helical fraction 



e = if:„A„... = |f . (7) 

71—1 

where the partition sum can be evaluated to give 

Z{N)^l + as^+'~{N + l)s+^^ . (8) 

(s- 1)^ 

2.1.5 The generalized zipper model 

The only difference is topological; then macromolecule consists of any num- 
ber of helical and coil regions which may alternate freely. One associates 
the following weights: 

• 1 if coil comes after helix or coil; 

• s if helix comes after helix; 

• as if helix comes after coil (nucleation). 



The model thus implements the ideas presented in section 2.1.2 without 
imposing any further constraints. The state of the residue at site i can be 
described by a 2-vector Vi, and the partition function is given by 

Zn ^ ^ < i'i\T\iy2 >< iy2\T\v3 > ... < vn-i\T\vn > 

{iyi}...{iyjv} 

J2 <>^i\T''Wn> (9) 

{'^i},{'^n} 

where the matrix elements of T express the Boltzmann factors specified 
above, i.e. 

s 1 
as 1 



(10) 



To evaluate the partition sum, I apply periodic boundary conditions (con- 
venient, not a must, and certainly wrong for small chains) and obtain 



Zn = TrT^ = + Af (11) 



where the eigenvalues are given by 



Ao,i = i[l + s±A] (12) 



A = V(l - s)2 + Aas 
and, in the large N limit, Z is dominated by the largest eigenvalue, Aq. 
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Note that the partition function (not the T -matrix!) can be mapped 
onto the one of the ferromagnetic Ising model with exchange interaction J 
and magnetic field /i, with the identifications 

a ^ (14) 

^magnetic — ^helix — coil • (-^^) 

To obtain the helix fraction, note that if the probability of obtaining a 
helical segment of length k is given by (j>k{<j)ks'^ , where 4> is the coefficient 
of s*^ in the partition sum. This gives 

e^--yM<^)ks'' = --^ (16) 

fc=i 

One can now verify that as s ^ 1, A ^ 6^1/2; for a << 1 (cf. 

below), 

'de\ _ 1 AH* 



(17) 



The experimental situation]^ for long (A^ = 1500) chains is summarized in 
Fig. [l]. Fits can be obtained with AH* — —3.8kJ/mole (cf. calorimetric 
measurements AH — —3. 97k J/ mole) and a = 1.6 x lO^"'. One usually in- 
terprets 1 / -y/CT as number of residues cooperatively involved in the transition 
(cf. AON theory, Eq. ||). This interpretation also follows from the Ising 
model, where the inverse correlation length is given (in lattice constants) 

by 

l/e = Ai - Ao = (at s = 1) . (18) 

The following simple conformational argument provides an independent 
estimate for a = g-AGi„it/T gASi„it_ initiation of the helix involves orga- 
nization of J (=3 or 4) residues, each one by 2 dihedral angles. Typically 
a dihedral angle can take 3 independent orientations in space. This gives a 
total of 3^"^ states, or an entropy loss ASmit — —2 Jin 3. (=-6.6 for J = 3, 
or -8.8 for J = 4). This compares favorably with ln(2 x lO"**) = —8.5. 

Similarly, one can relate the entropy loss involved in helix growth^ to the 
energy of the H-bond. At the transition, AS* = AH* /T^ = -1.85. This 
is roughly comparable to the estimate —2 In 3 ~ —2.20, obtained by con- 
sidering the 2 dihedral angles which must be organized to admit a residue 
into the helix. 



2.1.6 A useful shortcut 

It is possible to obtain the thermodynamics of the generalized zipper model 
without recourse to the transfer matrix formalism. I present this "hand- 
waving" S, because it will be useful for DNA loops (cf. below). 



6 



The fundamental "entity" of the macromolecule is a hehcal region of 
length n, followed by a coil region of length m. This "helix-coil" entity is 
characterized by a free energy 

g{n,m) = —Tina — nT Ins (19) 

where the two terms correspond to the contributions of helix nucleation 
and growth, respectively; it occurs with a probability 

Pn,m = exp{-[g(n,m) - go]/T} (20) 

where go = — Tlnz is the equilibrium free energy per site of the full macro- 
molecule, to be determined by the normalization condition 

oo 

Pn,m = l ■ (21) 

n,m— 1 

Both the n and the m- summations can be done trivially as long as s < z 
and z > 1. The condition ( pi] ) can then be written as 
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s z — s a 



(22) 



whose roots are identical to those of (|T^) obtained via the transfer matrix; 
the largest root is the one which satisfies the condition z > max(l, s) (cf. 
above) . 

2.2 Loop entropies and DNA denaturation 
2.2.1 Background 

Thermal DNA denaturation occurs when the two strands of the double 
helix separate upon heating. In real DNA the phenomenon of multistep 
melting is ubiquitous, reflecting the inhomogcncity of the molecule. "Ho- 
mogeneous", synthetic DNA, which consists of a few thousand identical 
base pairs has been studied experimentally |^ and shown to exhibit a very 
sharp transition; the qualitative shape of the coil fraction curve, as obtained 
by optical density or viscocity measurements, is similar to that of Fig. |l|; 
however, the observed temperature width is of the order of one degree. It 
is reasonable to speculate that in the thermodynamic limit the transition 
would be of the first order. Poland and Scheraga §(PS) proposed a sim- 
ple model of the thermodynamics involved, based on the ideas discussed 



in Section 2.1.2 , and the concept of loop entropy (cf. below). There is 
however a difference in the physical origin of the parameters involved; the 
helix growth probability s — exp{— e/T} now reflects the combined effect of 
both interactions which are significant in the bonded DNA state: the hydro- 
gen bonding responsible for binding base pairs and the stacking interaction 
between adjacent bases. 
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2.2.2 Loop entropies 



A denaturation loop of length m, i.e. a region of m sites where the double 
helix has locally melted, is characterized by the extra entropy it contributes 
This extra entropy has been calculated for lattice polymers and is of the 
form 

S'L(m) = am + 6 — clnm , (23) 

where a and b are constants and c depends on the dimensionality. In the 
case of Gaussian polymer chains, or random lattice walks, which ignore the 
effects of excluded volume, c = d/2. Taking account of excluded volume 
tends to increase the value of c. It will be seen below that this can have a 
decisive influence on the nature of the transition. 



2.2.3 The phase transition 

In an infinitely long DNA molecule, the fundamental entity is again a 
double- helical region of n sites (base pairs), followed by a denaturation 
loop of length m (2m bases); putting together the contributions from heli- 
cal and loop part (cf. Sections 2.2.1 and 2.2.2 ), I obtain the free energy of 
this entity as 



-l{n,m) 



-T\n<j + ne — mTlnu +Tc\nm 



(24) 



where u = exp(a/2) is another constant, and I have dropped the irrelevant 
constant b. It is now possible to derive the the thermodynamics exactly as 
in section 2.1.6| . Inserting ( p^ ) in the normalization condition ( pT| ) gives 



1 = aU 



(25) 



where s = exp(— e/T) is the only temperature dependent parameter, and 



t-(- 



U{x) = 

m— 1 

Near a; = 1 it is possible to approximate U(x) by the expressionjj 

n c-1 



(26) 



1 



1 



(27) 



where, for c> 1, f/o « f/i ~ C(c) 0. 

I now follow the thermodynamic behavior near the putative singularity 
by defining!^ an Sc = exp(— e/Tc) via 



u 

Sr. 



1 = aU{l) 



(28) 
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and subtracting (^8|) from (|2^) to obtain 



se - s + A + ^t/iA^-i = , (29) 

where 2 = tt^(l + A) and I have only kept lowest order terms in the small 
quantities A and s — Sc- It is now straightforward to use ( p9| ) and obtain 
A(s); the helix fraction is then given by 

91nz OA , , 

TT^CX— . 30 

a Ins OS 
Two cases can be distinguished: 

• 1 < c < 2. The linear term in A can be neglected in (p9|). The helix 
fraction is proportional to (Tc — T)'^'^~'^^l'^'^~^^ ^ i.e. it approaches zero 
continuously near the transition. In particular, if c = 3/2 (the value 
which corresponds to d = 3 and neglecting excluded volume effects), 
one obtains a second order transition. 

• c > 2. The linear term in A dominates, and the transition becomes 
first order, i.e. the helix fraction drops abruptly to zero at the tran- 
sition. 

The above analysis shows how crucial the value of c is in determining the 
nature of the transition. It has been long known[^ that excluded volume 
effects, as calculated within the framework of self-avoiding walks, can in- 
crease the value of c to 1.75 for loops embedded in three-dimensional space. 
Recent research Q suggests that c may be as high as 2.1. 



3 Hamiltonian approach to DNA denatura- 
tion 

3.1 The model 

A "minimal" Hamiltonian model of homogeneous DNA denaturation has 
been proposed by Peyrard and Bishop||l0|(PB). The model assumes two 
parallel, harmonic chains, with lattice constant Z, joined in the form of 
a ladder by anharmonic springs; the particular model proposes a Morse 
potential because of its analytical tractability, although any form with a 
repulsive core, a stable minimum and a flat top (e.g. Lennard- Jones) would 
be physically suitable. The emphasis is on modelling the unbinding of the 
two chains, not the helical aspect of the ordered state; this is done in the 
general spirit of the theory of critical phenomena, which has demonstrated 
that the "essentials" of the interactions completely determine the critical 
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behavior. The Hamiltonian 



Htot = — \U„ + V„ + UJq {Un - Un-l) + UJq {'"n - Vn-l) 

n 



(31) 



describes the motion of the two bound chains with coordinates {wn}, 
dots denote time derivatives; only the motion transverse to the chains is 
considered; bases have equal masses m and are connected by harmonic 
springs of equal strength, determined by the frequency Wq; the energy scale 
of the Morse potential 



V{x) = D{e- 



(32) 



is given by D and its spatial range by 1/a. 

Transformation to center-of-mass and relative coordinates, Yn = {un + 
Vn)/2,yn = Un — Vn, M = 2m,l/ jjL = 2/m decouples center-of-mass from 
relative motion, i.e. 

Htot = HoiY)+H{y) , 



where 



^ + \MuliYn-Yn..f 



Pn = MYn is the canonical momentum conjugate to Y^, and 



1 

^ + ^^^^o{yn - Vn-lf + V{yn) 



n 

where Pn = l^Vn is the canonical momentum conjugate to y. 



(33) 
(34) 

(35) 



3.2 Statistical Mechanics 

Hq is just the Hamiltonian of a harmonic chain with the total base pair mass 
2m per site. It gives an additive nonsingular contribution to all thermal 
properties. It will be neglected in what follows. The classical thermody- 
namics of H is described by the canonical partition function 



N 

Z= H dpndyne-f^" 

•' n=l 



(36) 



One can immediately perform the Gaussian integrals over momentum space 
and obtain 

Z = ZkZp, (37) 
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where each integration in the kinetic part contributes a {2t:^/ PY^'^ factor 
to the partition function, i.e. 

- (271^^7/3)^/2 . (38) 

The nontrivial part is 

Zp = j {^\\_dy,}^ T{yi,y2)---T{yN-i,yN)T(iiN,yN+i) , (39) 
where 



T{x,y) ^ 



(40) 



3.3 Transfer integral: the formalism 
3.3.1 Definitions and Notation 

Consider the eigenvalue problem defined by the asymmetric kernel T (the 
kernel can be easily symmetrized but need not be so; in fact, working with 
the asymmetric kernel is technically advantageous in examining the validity 
of some approximations, cf. below): 

dyT{x,y)'^^{y) = A,<i>^(x) (41) 

dyT{y,x)^iiy) = A,$^(a;) , (42) 

where left and right eigenstates have been assumed to be normalized; note 
that the normalization integral is J dx^^ {x)^^{x). Orthogonality 

dx^^{x)^^,{x)^S,,, (43) 

^ — OO 

and completeness 

Y,^i{x)^^'{y)^S{x-y) (44) 

relationships are assumed to hold. Note that this is not obvious for the 
class of potentials of interest here. This is a point which will be further 
taken in section 3.6, I will further use the notation 



K = e-^'^ (45) 
(sensible as long as the eigenvalues are nonnegative). 
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3.3.2 The partition function 



The integrand of (39), as written down has a problem: it includes a refer- 
ence to the displacement Un+i of the N + 1st particle, which has not yet 
been defined. For a large system, this is best remedied by means of periodic 
boundary conditions (PBC), i.e. by demanding that un+i — Vi- Alterna- 
tively, the integration may be extended to one more variable, dyisi+i, with 
the simultaneous introduction of a factor 5{yM+i — yi) to take care of PBC. 
This however is the same as the sum in the left-hand-side of (p4[). I then 
obtain 

Zp = y2 dyi - ■ ■ dyN+i'^Uyi)T{yi,y2) ■ ■ ■T{yN,yN+i)'^^{yN+i) ■ 

(46) 

The braces make clear that I can perform the integral over dj/Af+i and ob- 
tain a factor Ai^<i>^(yjv-(-i)j using the defining property of right-hand eigen- 
functions. The process can be repeated N times, each time giving a further 
factor Aj^ and a right eigenfunction with an argument whose index is smaller 
by one. At the end, I am left with 

Zp = Y.f dy^^^^{y,)K^^^{y^) = ^ A^^ . (47) 

In the thermodynamic limit, Zp is dominated by the largest eigenvalue Aq 
or, equivalently, the lowest eo: 

lim ^lnZp = In Ao = -/3eo (48) 

N—^oo i\ 

3.3.3 The order parameter 

<yi> = I dyi---dyNT{yi,y2)---T{yi_i,yi)yi 



Zp 

T{yt,yi+i) ■ ■■T{yN,yN+i) 



J '^Vi ■ ■ ■ dyN+i^y[yi) T{yi,y2) ■ ■ ■T{yi.i,yi) yi 



i-l 



T{y,, 2/,+i) • • • Tiy^, y^+i) ^^{yN+i) , (49) 



after insertion of a complete set of states (cf. above); the braces denote the 
number of times I can perform an integration and obtain, respectively, a 
right eigenfunction with an argument smaller by one, or a left eigenfunction 
with an argument larger by one, as well as a factor Aj,. The remaining 
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integral must be performed explicitly: 

« Moo (50) 

where the second line is exact in the thermodynamic limit, and I have used 
the abbreviation 

/oo 
dy^^{y)y^^iy) . (51) 
-oo 

3.3.4 Correlations 

With i < j, 

<VtVj> = J dyi---dyNT{yi,y2)---T{yi_i,yi)yiT{y„yi+i) 

■ ■ ■ T{yj-i,yj)yjT{yj,yj+i) ■ ■ ■ T{yN, y^+i) 
= -^^jd-Vi--- dyjK~^'^^{yi)yiT{yi, yi+i) 

•••T(y,_i,y,-)2/,-Ar^'+^*?fe) (52) 

where the straightforward integrations, i.e the first i—1 and the last N—j+1 
have aheady been performed (cf. above). In order to perform the remaining 
integrations, 1 insert two more factors of 1, after yi and before yj, i.e. 
integrals / 6{yi — iji) and / 5{yj — yj), respectively; exploiting the presence 
of the 5 functions, I may substitute the variables j/j and yj by y^ and yj 
respectively. This translates to two more sums over complete sets of states 
and another j — i integrals which can now be performed: 

< ViVj > = ^ E / dy^dyJdy^ ■ ■ ■ dyjAl-^^^{yi)y^^^^{y,)<i>l:iyi) 
T{yuyi+i) ■ • •T(j/,_i,t/,)$f (y,)$^(y,)y,AW$^(^.) 
= ^ E ^^^'-'^r j dyidy3^i{y»^{yi) 



\L r- \-. 



<5^,P*p (%■)%*.(%) 

= ^E^^^'"'^riM.^i' • (53) 



Zp 



In the thermodynamic limit the u = term dominates; the resulting factor 
cancels against the denominator and leaves 

< yiyi+r >= E I^omI' e-^^^"-^")'- (54) 
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where I have used Dirac shorthand for the matrix element and set j = i + r. 
The first term (/x = 0) in the above sum corresponds to < y >^ and should 
properly be subtracted from both sides; This leaves 

< 5yi5yi+r > = < ViVi+r > - < Vt >< Vt+r > (55) 

where now the ground state is excluded from the summation. The above 
result identifies the correlation length ^, i.e the typical length over which 
the decay of correlations takes place, as 

I P(ei-eo) 

where the subscript 1 stands for the first excited state (dominant exponen- 
tial in the limit of large r) . 



3.4 TI results: Gradient-expansion approximation 

Suppose that the displacement field does not change appreciably over a lat- 
tice constant. This is certainly reasonable at low temperatures. Note that 
this does not exclude large displacements per se. Nonlinearity is explicitly 
allowed, but the displacement field must be smooth. The assumption is 
certainly reasonable at low temperatures. 

I set y = X + z, <i>^ (p and rewrite ( ^l|) as 

= dz e-i'^M-o^-^ \4>,{x) + z(p'^{x) + \zm{x)\ 

J —OO 

27r 

where higher terms in the gradient expansion have been neglected and the 
Gaussian integrals have been performed; this is meaningful as long as the 
width of the Gaussians is smaller than the range of the Morse potential, 
i.e. 

Ptiujl/a^ > 1 . (58) 

The factor in front of the r.h.s. of ( |57| ) can be absorbed in the eigenvalue by 
defining €,y + l/{2(3)] ln[27r/(/3^a;o)]. Now, for many practical purposes, 
when it comes to calculating matrix elements, the relevant magnitude of 
e — V{x) is D, the depth of the Morse well (or some other characteristic 
energy in the case of another potential). The key to this statement is that 
one does not need to consider large negative values of x, where V{x) is huge, 
because at such x, both the exact eigenfunction $ and its approximation 



(t>v{x) 



,£{x) 



(57) 
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4> can be expected to be negligible. If then I3D < 1^ it is reasonable to 
expand the exponential and keep only the first term. Dividing both sides 
by /3, I obtain a Schrodinger - like equation, 

^ 4>l{x) + [V{x)-i,]Mx)^Q . (59) 



Before continuing the discussion of (|59|) and its properties, I pick up the 
bits and pieces (cf (^), (^, (|4^ ) ) of the thermodynamic free energy (per 
site) 

where f — eq. The first term in ( |60| ) is the free energy of the small os- 
cillations (transverse phonons in this context). It is a term smooth in 
temperature (constant specific heat!) and therefore irrelevant to any phase 
transition. Any nontrivial physics is hidden in the second term, which is 
identical with the the smallest eigenvalue of (|59). 
A couple of comments are in order. First, 



(59|) would be a literal (i.e. 
quantum-mechanical) Schrodinger equation, if I substituted l/(/3a;o) by h. 
I will come back to that point. Second, I can get a dimensionless potential 
(and eigenvalue) by dividing both sides of (|5^) by D. In other words, the 
relevant dimensionless parameter is 




(61) 



^ n + 1/ 2]' 



D (quantum mechanics) 

D (statistical mechanics). 
In terms of S, the bound state spectrum of (|59| ) is given ||l^ by 

fe - ' 

n ^ 0,l,...,mi((5- 1/2) . (62) 

There is at least one bound state if 6 > 1/2. For 1> S > 1/2 there is exactly 
one bound state. And if S becomes equal to, or smaller than 1/2, there is 
no bound state at all. The value 6c — 1/2 is "critical". In quantum 
mechanical language, if a particle has a mass which is lighter than a critical 
mass [ic = h'^a'^/{8D), it cannot be confined in the Morse well. Quantum 
fluctuations will drive it outPl In the context of statistical mechanics, 5c 



^Note that, in connection with fesh, this defines a temperature window D < kgT < 
/iaJn/a'^ for the vaUdity of the overall approximation scheme. 

^This is a general property of asymmetric one-dimensional wells; symmetric wells will 
support a particle in a bound state, no matter how low its mass. 
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corresponds, via (pi]), to a critical temperature Tc = 2{ujQ/a)^2fiD. The 
free energy is given by 



I 
D 



T>Tr 



(63) 

T " 



1 - ( 1 - ^ I T<T, 



where in the upper Unc I have made use of the fact that the bottom of 
the continuum part of the spectrum is at e = The free energy / is 
non-analytic at T — Tc, where its second derivative is discontinuous (i.e. 
there is a jump in the specific heat). This corresponds to a second order 
transition, according to the Ehrenfest classification scheme^. 

In order to gain some further insight into the physics involved^ it is 
useful to examine the average displacement (|5^), determined by the ground- 
state (GS) eigenfunction 

0o(x)=e-^/2C'-V2 (64) 

where ^ = 25e~°'^. It is straightforward to see that, as T approaches Tc 
from below, the eigenfunction extends towards larger and larger positive 
values of x: 

M^) « e"^^ (65) 

where 

is a (transverse) characteristic length which measures the spatial extent of 
the GS eigenfunction. As a consequence, we can estimate that < y >, 
which is dominated by the large values of the argument, will also behave as 

<y>^{S-6c)-'^ (^l-^^ ' . (67) 

As the critical temperature is approached from below, particles cease to 
be confined to the minimum of the Morse well. They perform larger and 
larger excursions to the flatter part of the potential. At Tc the transition is 

■^Note that the term "second order" is meant literally in this case, not just as a 
metaphor for the absence of a latent heat (for which the term "continuous transition" 
would be appropriate). 

* The mathematical analogy between the behavior of the spectral gap which occurs 
in a point (d = 0) system and the singularity in the free energy of a classical chain 
(d = 1) is an example of a deeper analogy which relates quantum to thermal fluctuations; 
the formal correspondence 7i <-> 1/ {ffuJo ) manifests a far-reaching analogy between d- 
dimensional quantum mechanics and {d+ l)-dimensional classical statistical mechanics. 
The analogy is most fruitful at d = 1, because of the interplay and the richness of exact 
available results which based either in the transfer-matrix approach of 2-dimensional 
classical statistics or on the Bethe-Ansatz developed for 1-d quantum spin systems. 
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complete; the average transverse displacement is infinite. Particles move, on 
the average, on the flat top of the Morse potential. Unwinding ( "melting" ) 
of the DNA has occurred. 

In the language of critical phenomena < y > is the order parameter. 
In the "usual" phase transitions, one goes from an ordered to a disor- 
dered phase. The order parameter m vanishes at the transition point, i.e 
m DC (Te — T)^ with a positive critical exponent (3 (not to be confused with 
the inverse temperature: standard notation of critical phenomena!). DNA 
melting is really an instability!^ - rather than an "order-disorder" transi- 
tion. It is therefore not surprising that the corresponding critical exponent 
/? extracted from ( |67|) is negative (-1). 

Experimental data on DNA denaturation do not deliver <?/> di- 
rectly. The "experimental order parameter" is the helical fraction, i.e. the 
probability that a given base pair is still bound; technically one uses an 
(instrumentation-dependent) cutoff yo and measures P{y > yo,T). For the 
model presented here, this function approaches zero smoothly (linearly) as 
T Tc, independently of the choice of t/o- 

Eq. (|5^) states that the correlation length is also contoUed by the gap 
in the eigenvalue spectrum; as the transition is approached, 

which identifies a critical exponent = 2 for the divergence of the corre- 
lation length. The picture of thermal denaturation which emerges is one 
of ordered regions, where helical structure persists; these regions are in- 
terrupted by droplets of the high-temperature phase, i.e. "denaturation 
bubbles" of typical size ^. 



3.5 A first order transition? 

It is possible to generalize the theory in order to take account of the fact 
that the stacking energy is a property of successive base pairs, rather than 
individual bases. A practical way of doing this is to substitute the second 
term in the Hamiltonian ( p5[ ) by 

i^cjg [1 + .g(y„ + ?/„_i)] (y„ - y„_i)^ , (69) 

where Q 

gix) = e-"^ . (70) 



The effect of Eqs. (|69|)-(70) is to interpolate between the original value 
of the elastic coupling if either (or both) of the two base pairs n, n — I 
is unbound (in which case ?;„ — > oo or yn~i — > oo), and twice that value 
if both are bound; in the latter case, typically, j/„ w 0; the much higher 
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values of the stacking energy, which (70) in principle allows, are statistically 



irrelevant due to the repulsive core of the Morse potential. Within the 
gradient expansion approximation, it can be shown that the main effect 
of the nonlinear stacking energy on the thermodynamics is to generate an 
effective, on-site, "thermally activated" barrier 

U{y) = ^ln{l + e-^-y) . (71) 

which appears in Eq. J59 and acts in addition to the Morse potential. 

It has been shown |p^p^, p7[ | that the character of the transition changes 
dramatically as the value of the stacking parameter ratio a/ a decreases 
(corresponding to a longer range in the effective potential). Although the 
transition remains asymptotically second order, the limiting asymptotic 
behavior becomes relevant only within an exponentially small range of the 
temperature difference — T. For all practical purposes, the transition 
is first order, with a finite melting entropy A5 — A^D /Tc, where is a 



numerical constant of order unity 1 15 



It should be noted that the interpolation ( |70D is not unic^ue; an inter- 
polation function of the type g{x) — \, x < xq, g{x) = (xo/x)"^, x > xq, 
leads - depending on the other parameters - to a rich variety of critical be- 
havior, ranging from a first-order transition to continuously varying critical 
exponents. 

With the above modification (|69|), it has become possible [T^ to de- 
scribe, at least in principle, the series of multistep melting observed in real, 
heterogeneous DNA. 

3.6 TI beyond the gradient expansion 

It was stated in Section that the TI formalism rests on the assumption 



that the integral equations (^l|) and ([42| ) - have a complete, orthonormal 
set of eigenfunctions. Within the gradient approximation approach this was 
demonstrated by construction - since the integral equation was reduced to a 
Schrodinger-like equation. In many cases however, the gradient expansion 
is not valid at all. It is therefore necessary to develop alternative, mostly 
numerical methods for computing TI thermodynamics. For such applica- 
tions it is expedient to consider the symmetrized, dimensionless version of 
the kernel (^), i.e. 

T,{x, y) = e-\^y-f/R+v{-)+v{y)]/{2T) (73) 
and the associated integral equation 

dyTs{x,y)^{y) ^ K4>{x) , (73) 
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where R = Da^ /{^ujq), T = 1/{PD) and the Morse potential is now dimcn- 
sionless, V{x) = (1 — e^"^)^, as are the displacement variables x,y. 

Due to the flat top of the Morse potential, the kernel ( [72|) is not of the 
Hilbert-Schmidt type|Q; therefore the integral equation (|73|) is singular 
and it can not be a priori stated that it possesses a complete orthonormal 
set of eigenstates; in other words, the prerequisites for directly applying the 
TI method are not strictly met. In the rest of this section I will outline a 
mathematically consistent procedure of examining the spectral gap of ( fz^ ) , 
based on finite-size scaling concepts 

Due to the presence of the Gaussian factors in the kernel, it is possible 
to approximate the integral in the left-hand-side of (|73| ) by using a Gauss- 
Hermite grid of size A^, i.e. 

N 

dye-y'f{y)^Y.'"rnf{yra) (74) 

• m=l 

where the positions {j/m} and weights {wm} are given by the appropriate 
Gauss- Hermite quadratures routine. The largest y^r ~ {IN -1-1)^/^ = L 
can be used as estimate of the transverse "size of the system" employed at 
any given discretization. I emphasize transverse because the length of the 
chain is infinite, i.e. the thermodynamic limit has already been taken. 

I use "rescaled" variables, i.e. y — py, p — (2RTy/^, divide both 
sides of (^ by pV^, and use the approximation (|7j). The result is an 
approximation of ( |73| ) by the matrix eigenvalue equation 

N 

Y,D,,A'; = KA'i (75) 



where 



l^mWf^ 1/2 ^y'yO ^-(^yi-y^f/2^-[V{pr) + V{py^)]/i2T) ^-yg^ 



and A, = A,/(27ri?T)i/2 = e"^-/^. 

It is now possible to solve numerically the real, symmetric matrix eigen- 



value problem (76) for a range of temperatures and a sequence of increas- 
ingly fine grids. Results for the difference between the two lowest eigenval- 
ues are shown in Fig. || for R = 10.1 . For any given size L, the gap has 
a minimum Ae,„(i) at a certain temperature Tm{L). Fig. ^demonstrates 
that (i) the value of the gap approaches zero quadratically as L — > cx3 (with 
an accuracy of 10^^), and (ii) the sequence of T„i{L)'s also approaches a 
limiting value Tc — 1.2275 quadratically. 

It is natural to identify the limiting temperature Tc, where the spectral 
gap of the limiting, infinite-dimensional matrix eigenvalue equation ((Tq ) 
vanishes, as the transition temperature of the original TI equation (|73|). 



19 



0.08 



0.00 



0.06 



0.04 



0.02 




0.9 1.0 1.1 1.2 1.3 
T 

Figure 2: The gap between the two lowest eigenvalues of the matrix eigen- 
value problem (^5|), for a variety of N values. For a given iV, the gap has 
a minimum at a certain temperature Tm- 

Further application of finite-size scaling methods demonstrates that the 
various critical exponents coincide with those obtained within the gradient 
expansion method. 
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Appendix: Phase transitions in one-dimensional 

systems 

I briefly discuss why the general prohibitions on phase transitions in one 
dimension are inapplicable to both the PS and the PB models of DNA 
denaturation. 

Van Hove's theorempT| states that no phase transitions occur in 1-d 
particle systems with short-range pair interactions. The PB model has on- 
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Figure 3: The magnitude of the gap minimum (circles, right y-axis scale) 
approaches zero as the system size goes to infinity. The sequence of the 
temperatures corresponding to the gap minima, Tm{L) (diamonds, left y- 
axis scale), can be used to provide an estimate of the critical point Tc- 



site potential - i.e. the theorem is not applicable. It is however worth noting 
that similar mathematical proofs, have been given for systems with periodic 
on-site potentials |2^. Such proofs however seek to prove analyticity of the 
eigenvalue spectrum and hence absence of a phase transition; as such, they 
tend to exclude potentials which give rise to singular TI equations. 

The PS model is not a Hamiltonian model and therefore van Hove's 
theorem is again not applicable. 

Landau's theoremps] is significantly stronger. It states that "macro- 
scopic phase coexistence cannot occur at finite temperatures in one dimen- 
sional systems" . It is less obvious why it should not apply. I therefore 
outline the proof. Consider a system with N sites, which may exist in ei- 
ther phase A or phase B. Let 6 be the fraction of phase A; furthermore, 
let there be m << iV contacts between the phases, each of energy e. These 
can be steplike (Ising) or continuous domain walls. The free energy of the 
configuration is given by 

F = NOJa + ^(1 - e)fB + Fdw (A.77) 
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where 

Fdw = me — kBTSDwini, N) 
and the (dimensionless) entropy is given by 



Sowim, N) — In 



m!(iV - m)! 



dn 



Ne 
m 



(A.78) 



(A.79) 



Minimization of the total free energy with respect to m yields a macroscopic 
average number (ie. a finite density) of domain walls 



(A.80) 



The system breaks up into m regions of finite size e'^^'^'^'^'^K Macroscopic 
phase separation can only occur at zero temperature (as the domain size 
goes to infinity). 

Landau's argument covers a wide range of systems, e.g. double-well 
on-site potentials (Ising universality class), or periodic on-site potentials. 
It does not cover the PB case, because the DW has infinite energy It 
does not apply to the PS case becau se the loop entropy is not proportional 
to the size of the loop and therefore ( A. 77 ) does not hold. On the contrary, 
the theorem is applicable to the generalized zipper model, as its authors 
had correctly noted [Q. 
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